###################################################
## SonPark_replication_SI.R
## Replication R code for SI
## Date: "Thu Aug 11 10:45:26 2022"
###################################################
rm(list=ls())
start_time <- Sys.time()
###################################################
## Load packages
###################################################
library(ggplot2)
library(tibble)
library(gridExtra)
library(readxl)
library(dplyr)
library(RColorBrewer)
library(ggthemes)
library(kableExtra)
library(broom)
library(jtools)
library(car)
library(ggthemes)
library(stargazer)
library(xtable)
library(ggh4x)
library(gridExtra)
library(visreg)
library(dplyr)
library(arsenal) 
library(table1)
library(Rmisc)
library(tidyr)
library(sjPlot)
library(table1)

###################################################
## Load data
###################################################
load("SonPark_Replication.RData")

## Treatment level information
Treatment.level.2019 <- c("Conditional Military Punishment",
                     "Reduced Nuclear Threat",
                     "Economic Sanction",
                     "Resolution on ROK's NPT Violation",
                     "Nuclear Technology Sanction",
                     "Elite Opposition to Proliferation",
                     "Public Opposition to Proliferation")
df2019$Treatment <- factor(df2019$Treatment, levels=Treatment.level.2019)
df2019$Treatment <- relevel(df2019$Treatment, ref = "Public Opposition to Proliferation")

Treatment.level.2022 = c("Conditional Military Punishment",
                           "Reduced Nuclear Threat",
                           "Economic Sanction",
                           "Resolution on ROK's NPT Violation",
                           "Nuclear Technology Sanction",
                           "Elite Opposition to Proliferation",
                           "Public Opposition to Proliferation",
                           "Control")

df2022$Treatment <- factor(df2022$Treatment, levels=Treatment.level.2022)
df2022$Treatment <- relevel(df2022$Treatment, ref = "Control")

###################################################
## SI Table 1 
###################################################
#### Randomization Check for 2019 Experiment 
fit1 <- glm(Treatment=="Conditional Military Punishment" ~ age + male + ideology+political.knowledge+income+education, data=df2019)
fit2 <- glm(Treatment=="Reduced Nuclear Threat" ~ age + male + ideology+political.knowledge+income+education, data=df2019)
fit3 <- glm(Treatment=="Economic Sanction" ~ age + male + ideology+political.knowledge+income+education, data=df2019)
fit4 <- glm(Treatment=="Resolution on ROK's NPT Violation" ~ age + male + ideology+political.knowledge+income+education, data=df2019)
fit5 <- glm(Treatment=="Nuclear Technology Sanction" ~ age + male + ideology+political.knowledge+income+education, data=df2019)
fit6 <- glm(Treatment=="Elite Opposition to Proliferation" ~ age + male + ideology+political.knowledge+income+education, data=df2019)
fit7 <- glm(Treatment=="Public Opposition to Proliferation" ~ age + male + ideology+political.knowledge+income+education, data=df2019)
stargazer(fit1, fit2, fit3, fit4, fit5, fit6, fit7)


###################################################
## SI Table 2 
###################################################
#### Randomization Check for 2022 Experiment 
fit8 <- glm(Treatment=="Control" ~ age + male + ideology+political.knowledge+income+education, data=df2022)
fit9 <- glm(Treatment=="Conditional Military Punishment" ~ age + male + ideology+political.knowledge+income+education, data=df2022)
fit10 <- glm(Treatment=="Reduced Nuclear Threat" ~ age + male + ideology+political.knowledge+income+education, data=df2022)
fit11 <- glm(Treatment=="Economic Sanction" ~ age + male + ideology+political.knowledge+income+education, data=df2022)
fit12 <- glm(Treatment=="Resolution on ROK's NPT Violation" ~ age + male + ideology+political.knowledge+income+education, data=df2022)
fit13 <- glm(Treatment=="Nuclear Technology Sanction" ~ age + male + ideology+political.knowledge+income+education, data=df2022)
fit14 <- glm(Treatment=="Elite Opposition to Proliferation" ~ age + male + ideology+political.knowledge+income+education, data=df2022)
fit15 <- glm(Treatment=="Public Opposition to Proliferation" ~ age + male + ideology+political.knowledge+income+education, data=df2022)
stargazer(fit8, fit9, fit10, fit11, fit12, fit13, fit14, fit15)


###################################################
## SI Table 6
###################################################
#### (1) logistic regression on opinion change ####
fit1 <- glm(opinion.change ~ Treatment, data=df2019, family = "binomial") ;summary(fit1)
fit2 <- glm(opinion.change ~ Treatment + age + male + ideology  + political.knowledge +
              income + education, data=df2019, family = "binomial") ;summary(fit2)

#### (2) logistic regression on behavior change ####
fit3 <- glm(behavior.change ~ Treatment, data=df2019, family = "binomial") ;summary(fit3)
fit4 <- glm(behavior.change ~ Treatment+ age + male + ideology  + political.knowledge +
              income + education, data=df2019, family = "binomial") ;summary(fit4)
fit5 <- glm(behavior.change ~ Treatment + relevel(choice.at.treat, ref="No Change") + age+ male + ideology  + political.knowledge +
              income + education, data=df2019, family = "binomial") ;summary(fit5)

#### (3) logistic regression on attitude change ####
fit6 <- glm(attitude.change ~ Treatment, data=df2019, family = "binomial") ;summary(fit6)
fit7 <- glm(attitude.change ~ Treatment + age + male + ideology  + political.knowledge +
              income + education, data=df2019, family = "binomial") ;summary(fit7)

stargazer(fit1, fit2, fit3, fit4, fit5, fit6, fit7)


###################################################
## SI Table 7
###################################################
#### (1) logistic regression on opinion change ####
fit1 <- glm(opinion.change ~ Treatment, data=df2022, family = "binomial") ;summary(fit1)
fit2 <- glm(opinion.change ~ Treatment + age + male + ideology  + political.knowledge +
                income + education, data=df2022, family = "binomial") ;summary(fit2)

#### (2) logistic regression on behavior change ####
fit3 <- glm(behavior.change ~ Treatment, data=df2022, family = "binomial") ;summary(fit3)
fit4 <- glm(behavior.change ~ Treatment + age + male + ideology  + political.knowledge +
                income + education, data=df2022, family = "binomial") ;summary(fit4)

#### (3) logistic regression on attitude change ####
fit6 <- glm(attitude.change ~ Treatment, data=df2022, family = "binomial") ;summary(fit6)
fit7 <- glm(attitude.change ~ Treatment + age + male + ideology  + political.knowledge +
                income + education, data=df2022, family = "binomial") ;summary(fit7)

stargazer(fit1, fit2, fit3, fit4, fit6, fit7)

###################################################
## SI Table 8
###################################################
##############################################################################
## changing the reference to other types of nonproliferation information 
##############################################################################
fit.ct <- list()
for(j in 1:7){
  df2019$Treatment <- relevel(df2019$Treatment,
                                ref = Treatment.level.2019[j])
  fit.ct[[j]] <- glm(attitude.change ~ Treatment +  ideology  + political.knowledge +
                       income + education + male + age, data=df2019, family = "binomial") 
}
## https://stats.oarc.ucla.edu/r/modules/coding-for-categorical-variables-in-regression-models/
stargazer(fit.ct[[1]], fit.ct[[7]], fit.ct[[2]], fit.ct[[3]],  fit.ct[[4]], fit.ct[[5]], fit.ct[[6]], no.space=TRUE)
## Reorder the rows so that it matches the printed version
## It may not work in some machines. Then, print them separately.
## ex) stargazer(fit.ct[[1]], fit.ct[[2]], fit.ct[[3]],  fit.ct[[4]], fit.ct[[5]], no.space=TRUE)
## stargazer(fit.ct[[6]], fit.ct[[7]], no.space=TRUE) ## fit5, fit6, fit5, fit6, 


###################################################
## SI Table Figure 5 
###################################################
require(quanteda)
library(readtext)
set.seed(10000)
df.back <- readtext("backfire.txt")
#### Use "Papago" for translating KOR to ENG ####
#### https://papago.naver.com/ ####

corp <- corpus(readtext("~/Dropbox/Nuclear/Data/backfire.txt")) 
corp_line <- corpus_segment(corp, "\n\n",  extract_pattern = FALSE, pattern_position = "after")

corp_line <- tolower(corp_line)
corp_line <- gsub("self-defense", "selfdefense", corp_line)
corp_line <- gsub("north korea", "northkorea", corp_line)
corp_line <- gsub("united states", "us", corp_line)
corp_line <- gsub("nuclear weapon", "nuclearweapon", corp_line)


# remove all digits
(corp_supply <- stringi::stri_replace_all_regex(corp_line, "\\d", ""))
## remove all punctuation and symbols
(corp_supply <- stringi::stri_replace_all_regex(corp_line, "[\\p{p}\\p{S}]", ""))

toks <- tokens(corp_supply)  %>%
  ## remove_punct = TRUE, 
  ## remove_symbols = TRUE, 
  ## remove_number = TRUE) %>%
  tokens_select(min_nchar = 2)
print(toks)

## compound
toks <- tokens_remove(toks, pattern = stopwords("en")) ## %>%
## tokens_ngrams(n = 2:3)
## toks_ngram

dfmat <- dfm(toks)
(dfmat.stem <- dfm_wordstem(dfmat))
require(quanteda.textstats)
require(quanteda.textplots)
require(quanteda.corpora)

g1 <- dfmat %>% 
  textstat_frequency(n = 10) %>% 
  ggplot(aes(x = reorder(feature, docfreq), 
             y = docfreq)) +
  geom_point() +
  coord_flip() +
  labs(x = NULL, y = "Frequency") +
  theme_minimal()

ggsave(g1,filename=paste0("SI_Figure5.pdf"),
       device=cairo_pdf,  
       width=6, height=3, family="sans")


###################################################
## SI Table 9
###################################################
tab11.2019 <- with(df2019, table(choice.at.treat, behavior.change))
df.tab11.2019 <- data.frame("No" = tab11.2019[,1],
                       "Yes" = tab11.2019[,2],
                       "Yes/Row total" = tab11.2019[,2]/apply(tab11.2019, 1, sum))
xtable(df.tab11.2019)


end_time <- Sys.time()
print(end_time - start_time)
